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Abstract 

Bimolecular diffusion coefficients are important parameters used by 
atmospheric models to calculate altitude profiles of minor constituents in 
an atmosphere. Unfortunately, laboratory measurements of these coef¬ 
ficients were never conducted at temperature conditions relevant to the 
atmosphere of Titan. Here we conduct a detailed uncertainty analysis 
of the bimolecular diffusion coefficient parameters as applied to Titan’s 
upper atmosphere to provide a better understanding of the impact of 
uncertainty for this parameter on models. Because temperature and pres¬ 
sure conditions are much lower than the laboratory conditions in which 
bimolecular diffusion parameters were measured, we apply a Bayesian 
framework, a problem-agnostic framework, to determine parameter esti¬ 
mates and associated uncertainties. We solve the Bayesian calibration 
problem using the open-source QUESO library which also performs a prop¬ 
agation of uncertainties in the calibrated parameters to temperature and 
pressure conditions observed in Titan’s upper atmosphere. Our results 
show that, after propagating uncertainty through the Massman model, 
the uncertainty in molecular diffusion is highly correlated to temperature 
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and we observe no noticeable correlation with pressure. We propagate 
the calibrated molecular diffusion estimate and associated uncertainty to 
obtain an estimate with uncertainty due to bimolecular diffusion for the 
methane molar fraction as a function of altitude. Results show that the 
uncertainty in methane abundance due to molecular diffusion is in general 
small compared to eddy diffusion and the chemical kinetics description. 
However, methane abundance is most sensitive to uncertainty in molecu¬ 
lar diffusion above 1200 km where the errors are nontrivial and could have 
important implications for scientific research based on diffusion models in 
this altitude range. 


1 Introduction 

Titan atmospheric models contain two major modules: a chemical kinetics mod¬ 
ule and a transport module. While there have been numerous efforts to quantify 
uncertainty in the chemical kinetics module [13, 14, 3, 5, 2, 9, 15, 23, 22, 11], the 
transport module is not as thoroughly evaluated. Furthermore, the transport 
module consists of two components: an eddy diffusion component and a molec¬ 
ular diffusion component. [4] assessed the question of uncertainty in chemical 
transport, but only focussed on the eddy diffusion profile effects, since it is the 
main contributor to uncertainty in the lower atmosphere, [34, Sec. 2.2.3] deter¬ 
mined by a Monte Carlo approach a best fit of eddy diffusion profile. However, 
molecular diffusion plays a critical role in determining altitude profiles in the 
upper atmosphere and has yet to be investigated cognizant of uncertainty and 
sensitivity. Parameters for molecular diffusion—molecular binary diffusion co¬ 
efficients for each pair of molecules—have not been revisited since their initial 
measurements made in 1973 [32]. 

[35] reviewed binary diffusion, describing three different sources for the bi¬ 
nary molecular diffusion coefficients, all of which are the result of a nonlinear 
least-squares fit to experimental measurements. These measurements were per¬ 
formed under temperatures ranging from 300 K to 700 K at 1 atm (see Fig. 1). 
In order to model bimolecular diffusion in Titan’s atmosphere, measurements 
must be extrapolated to conditions relevant to Titan—approximately 150 K 
and 10~ 6 atm. Here we evaluate the potential error involved in propagating to 
temperature and pressure relevant to Titan’s atmosphere by also propagating 
its uncertainty through this extrapolation. 

In this work we apply a probabilistic calibration framework based on Bayes’ 
theorem. The application of this Bayesian framework solves the problem of es¬ 
timating, with uncertainty, unknown molecular diffusion parameters in a model 
for Titan’s atmosphere. We review a brief history of the evolution of bimolecular 
diffusion models in section 2, and in section 3 we set up the Bayesian framework 
in which we estimate parameters for bimolecular diffusion. In section 4 we de¬ 
tail the measurements and measurement error of bimolecular diffusion available 
in the scientific literature and we couple these measurements with a Bayesian 
calibration procedure to estimate parameters for bimolecular diffusion. This 
estimate is then propagated through a bimolecular diffusion model to obtain an 
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estimate of bimolecular diffusion of the N 2 /CH 4 pair in physical conditions con¬ 
sistent with those typically found on Titan’s atmosphere. Using this, we then 
estimate the mixture diffusion coefficient of N 2 . We execute this entire process 
while taking into account all available experimental and a priori uncertainty. 
Section 5 concludes and summarises the paper. 


2 Molecular diffusion models 


Molecular diffusion in Titan’s atmosphere is modeled using the Wilke equation 
[33, 35], or the mixture-average rule, from the bimolecular diffusion coefficient: 
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where ntotai and n s are the total molar concentration and molar concentration of 
species s, respectively, and _D S is the bimolecular diffusion coefficient between 
species s and species j rn . The j m notation is used to take into account a possible 
approximation of the atmosphere’s composition where only the dominant species 
are considered, thus defining a “medium” in which the molecules are diffusing. 
In Titan’s atmosphere, this medium is typically composed of N 2 , CH 4 and 
sometimes H 2 . 

[8] derived a modified version of the mixture-averaged rule applicable for gas 
mixtures with more than two components. The purpose of this derivation was 
to address the rapidly changing composition of Titan’s upper atmosphere above 
the homopause, where binary diffusion is dominant. It required a rederivation of 
the transport term due to molecular diffusion, adapted to minor species diffusing 
through a bulk gas, 
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where H s is the scale height of species s, ar is the thermal coefficient, n to tai is 
the total density of the atmosphere and 
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where in s is the molecular mass of species s and m^ s is the mean molecular mass 
of the atmosphere without species s. This expression of molecular diffusion is 
suited for both minor and major species, with Eq. 3 converging towards Eq. 2 
for the minor species. The detailed analysis is given in [8, Section 3.3.1]. 

Three different models have been proposed in the literature for bimolecular 
diffusion between any pair of molecules, denoted as the pair of molecules s and 
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Figure 1: Measurements and fits. The uncertainty in the laboratory measure¬ 
ments is 4%. Titan’s atmospheric conditions are highlighted to emphasize the 
extrapolation required to apply these values to modeling Titan’s atmosphere. 


to. Note that a binary diffusion coefficient is specific to a pair of molecules and 
is agnostic to which species is diffusing and which species is the medium. First, 
[32] proposed a model depending only on temperature (Eq. 5), 

D a , m (T) = AT S , (5) 

where T is the temperature and A and s are the model parameters. Note that the 
temperature should be a reduced unitless temperature: it is more appropriate to 
write instead of T with T re f = 1 K. The second model, by [35], introduced 
a pressure dependence, 

D s , m (T,P) = —TP, (6) 

^-tot 

where ritot is the total molar concentration, A and (3 are the parameters. The 
pressure dependence is given by the mixture equation of state which, in the 
case of Titan’s atmosphere, is the ideal gas law, n to t = P/RT. Note that, 

just as above, it is more appropriate to write this with respect to a reference 

temperature. Finally, [20] also proposed a model with explicit temperature and 
pressure dependence, but included the bimolecular diffusion, D^°\ at 0°C and 
1 atm (T° and P° respectively) as a reference value, 

pO / j 1 \ 0 

D a , m (T,P) = DW—\^ J . (7) 

with and (3 being the parameters of the model. 
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Modeling Titan’s atmosphere requires a pressure extrapolation, it is there¬ 
fore more appropriate to consider either the model of [20] or [35]. Those models 
are equivalent through the state equation. From a modeler’s point-of-view, the 
model of [20] has the advantage to extrapolate a reference value of the diffusion 
coefficient, which holds therefore the potential to be informed independently. 

The rest of this work uses the model of [20]. 


3 Bayesian framework 

In many scientific applications, some quantities need to be determined from 
experimental observations that were subject to errors during the experiment. 
Since the quantity is an extrapolation of measurements, an exact value of this 
quantity cannot be known with certainty. It is therefore important to provide 
a measure of the associated uncertainty. We do this for the bimolecular dif¬ 
fusion model described above in a Bayesian setting. The Bayesian setting is 
advantageous because it provides a distribution of a quantity of interest, and 
this is achieved in two steps. The first step is calibration, where the model is 
fit to the data to obtain distributions for model parameters. The second step 
propagates these distributions to a quantity of interest. The Bayesian approach 
is not a statistical tool, but rather a method for describing an inverse problem. 
It differs from the simple Monte Carlo technique commonly used in planetary 
science applications. In the simple Monte Carlo technique one would use a 
Monte Carlo model as a tool to cast random experimental measurements from 
the measured mean values and errors, and then apply a Massman model fitting 
to each set of randomly determined measurements in order to obtain a probabil¬ 
ity distribution of D or /3. In the Bayesian approach, we use a Markov chain 
Monte Carlo (MCMC) algorithm (a statistical tool) to solve a Bayesian inverse 
problem. This approach is very similar to the simple Monte Carlo technique 
except that the MCMC does not cast random observations, as described above, 
but D and /3. These proposed values are then compared to the observations 
and the Metropolis-Hastings acceptance probability decides whether or not the 
proposed sample should be kept. The samples of and /3 at the end 

of the procedure are realizations from the probability distribution p(D^°\/3\y), 
which is desired output for our study. 

For the calibration step, we need to describe the model to calibrate and the 
data to calibrate against. We will calibrate the Massman model for bimolecular 
diffusion (Eq. 7) against direct noisy observations of D s m at various tempera¬ 
tures and pressures, 

Ujk D St m{Tj, Pk) T Vjki Vjk ^-A/"(0, CTjfc), (8) 

where j = 1,..., J, and k = 1,..., K. For the sake of simplicity, we write the 
observations as a column vector, y = (yn,..., yjx) T ■ Here Af{0, a 2 ) is notation 
for a Gaussian distribution with mean zero and variance a 2 . The a-y- will appear 
later in the expression for the likelihood (see Section 3.2). 
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The calibration step proceeds by finding the ‘best’ model parameters, 
and /3, given the observations y. That is, we seek the joint probability distribu¬ 
tion p((3, D^ly). Applying Bayes’s theorem yields, 

pQ3,D {0) \y\ oc p(y\/3^D < ' 0) )p{P, £> (0) ) . (9) 

posterior likelihood prior 

The prior distribution encodes a previously held state of knowledge, or an ex¬ 
ternal expert opinion, about what the parameters /? and D should look like. 
It is therefore given. The likelihood distribution is also known; given /? and 
£>(°) 

, one inserts them into the Massman model and evaluates Eq. 8. Therefore, 
since the right-hand side of Eq. 9 is known, the left-hand side is also known. 
The posterior distribution is the solution to the Bayesian calibration problem 
and there are many numerical methods to understand its properties. We choose 
to approximate the posterior distribution by statistical samples (see Section 
3.3). Although we have applied the Bayesian calibration framework to a spe¬ 
cific model, in general the framework is problem-agnostic and can be used in 
any scientific domain for parameter estimation and uncertainty quantification. 

The second step is to propagate the posterior distribution above, which we 
choose to approximate by samples, to the quantity of interest. For our purposes, 
the quantity of interest is the bimolecular diffusion D StTn (T, P) at some desired 
temperature T and pressure P. This is very easily done by evaluating 
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for each sample n = 1,..., N. This will yield another sequence of samples of 
bimolecular diffusions at the desired temperature and pressure. An estimate 
of the bimolecular diffusion at temperature T and pressure P can be obtained 
by computing the sample mean. The associated uncertainty of this estimate is 
obtained by computing the sample variance. 


3.1 Choice of prior 

The choice of a prior distribution is a very important part of the Bayesian 
framework, as it has a direct influence on the posterior’s description (see Eq. 9). 
In general, this choice can be very difficult. It is a statement about knowledge of 
the solution of the problem, which is unknown. However there are situations in 
which one can provide a prior distribution without knowing the solution exactly, 
but knowing some property of the solution. For example, an external observer 
can never know the exact speed of a car but it is assuredly non-negative and 
less than the speed of light. Although contrived, this example illustrates that 
common sense can be harnessed when expert guidance is absent. 

When a domain expert is present, one may attempt to acquire high-quality 
information which can be incorporated into a prior. Of course, the aim here is 
to do this as accurately as possible. This is still an active area of research for 
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which no gold standard exists [17]. The process of obtaining expert opinion for 
a particular problem is called elicitation. 

When one truly possesses no problem insight, the prior typically relies on 
the maximum entropy principle (MEP) [29, 10]. This principle states that, 
given a set of constraints, the distribution that encodes the least amount of 
information is one which maximizes the Shannon entropy [30]. For instance, 
if a mean and standard deviation of a parameter are given as constraints, the 
distribution that maximizes the Shannon entropy is a Gaussian. This method 
provides a distribution that attempts to maximize ignorance about the solution 
of the problem. Of course, when a domain expert is available one can still 
utilise the maximum entropy principle by setting the constraints through expert 
elicitation. 

The parameter is a diffusion coefficient and must be positive whilst 
/3 is an exponent and therefore can take any real value. With no knowledge 
other than bounds, the probability distribution that maximises ignorance in the 
Shannon entropy sense is a uniform distribution. Therefore, a priori, we choose 
£)(°) ~ t/[0,oo] and /? ~ [/[—oo,oo] so the prior distribution is a joint uniform 
on the /3) pair. 

3.2 Likelihood 

In Eq. 9 we must evaluate the likelihood of observing the data y given values of 
the parameters and /?. Since the experimental data was taken at a fixed 
set of temperatures T £ R J and pressures P € R^, we define the function Q to 
be, 


g(D<V,P) = D. itn { T,P) 

1 

TO ) 


p 


where division, multiplication and exponentiation are all done component-wise 
on the vectors T and P. We then derive the likelihood by recalling from Eq. 8 
that the difference of the model output and the observation is a Gaussian random 
variable. Therefore, the likelihood distribution function is given by a Gaussian 
PDF, 

p(y\f3, | exp i (s(/3, £><°>) - y) T E" 1 (s(/3, D™) - y)) . (11) 

The coefficient Z is a constant of proportionality and need not be computed. 

The matrix E is called the error covariance matrix and is equal to diag(txn,..., <Jjk), 
where crjk is as in Eq. 8. In general the expression of the likelihood function 
is determined by the knowledge available about the probability distribution of 
the experimental errors. Since most experimental results consist of an average 
of many measurements, and averages typically follow a Gaussian distribution, 
it is commonplace to express the likelihood as a Gaussian [16]. 
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3.3 Markov chain Monte Carlo 


The Bayesian framework described above is nothing more than a general prob¬ 
lem statement expressing the distribution of model parameters given noisy ob¬ 
servations in terms of a prior and a likelihood. To solve problems posed in 
this framework, numerical methods are used to understand the shape of the 
posterior, its mean, its variance, or to compute probabilities. We choose to 
use a method called Markov chain Monte Carlo to understand the posterior 
distribution on /3 and D 

Given a general probability distribution function p(x), Markov chain Monte 
Carlo (MCMC) methods produce a sequence of samples that, when plotted via a 
histogram, approximate p{x). In the limit of infinite samples, the approximate 
distribution converges to p. Most MCMC methods rely on the Metropolis- 
Hasting algorithm [21]. 

In our case p will be the posterior distribution function in the Bayesian 
framework (Eq. 9) and the MCMC will result in N samples, {/3 n , 
from the posterior distribution. One can compute the mean and variance of 
these samples to obtain an estimate with associated uncertainty. There are 
many variants of this algorithm each with their own benefits and drawbacks 
[6, 28, 1, 26, 27]. Initialising the Markov chain is also a challenging task, since 
a bad initial choice can be detrimental to sample quality [18, 24]. Typically, the 
initial condition is a guess provided by the user and therefore may lie in a low 
probability region with respect to the posterior distribution. In what follows, 
we maximise the posterior distribution to find a MAP (maximum a posteriori) 
point before starting the sampler. This ensures the Markov chain starts in 
stationarity. 

Samples in this work were generated by the delayed-rejection adaptive- 
Metropolis MCMC algorithm (DRAM) [12] implemented in the QUESO library 
[25]. QUESO is a C++ library for quantifying uncertainty in Bayesian calibration 
problems. It supports large-scale models that use many processors, and utilises 
multi-core architectures to provide high-quality samples from probability dis¬ 
tributions using Markov chain Monte Carlo. QUESO is free and open source 
software, available at http://libqueso.com. 


4 Results 

4.1 Data available 

The data available for the diffusivity of the (CH 4 ,N 2 ) couple are the measure¬ 
ments and fits of [32] and the fits of [20]. 

It is clear that the prior distribution influences the resulting posterior dis¬ 
tribution (see Eq. 9). We chose a minimally informed prior distribution, ob¬ 
tained by solving the maximum entropy principle. We assert that D 1°) is real 
and non-negative and (3 is real. This provides uninformative uniform priors 
ZA°) ~ U[0, oo] and /3 ~ U[— oo,oo]. 



T(K) 

Pch 4) n 2 (cm 2 s 4 ) 

Relative error (%) 

313.7 

0.242 

4 

314.9 

0.250 

4 

375.2 

0.353 

4 

474.7 

0.542 

4 

481.0 

0.539 

4 

573.5 

0.720 

4 

671.1 

0.919 

4 


Table 1: Data from [32] used for the calibration. 


[32] performed their measurements with pressure P(wa) = 1 atm. Using the 
ideal gas assumption, we have the following conversion equations, 


/?(Ma) — /3(Wa), 

D m _ , p (w.) ( T° \ 


( 12 ) 


where A( Wa ) and /3(wa) are the [32] model parameters, /?(Ma) and the [20] 
model parameters. 


4.2 Modeling conditions 

The temperature profile and medium species density conditions under which 
our calibration and uncertainty propagation have been performed were chosen 
to match those derived from data obtained during Cassini’s 40th flyby of Titan, 
commonly referred to as T40. The temperature profile, shown in Fig. 2, is based 
on fits of diffusion models to density measurements made by instruments on the 
Cassini orbiter and the Huygens lander as described in 

4.3 Calibration and convergence of calibrated quantities 

Fig. 3 shows 10 4 samples from the posterior distribution of D To assess how 
adequately a finite number of samples represent the posterior distribution, one 
typically investigates the convergence of the moments of the sampled variables 
with respect to sampler iteration. The moments we are interested in are the 
mean and variance of \ y and (3\. y (Fig. 4). 

Fig. 4 shows the mean (E) and variance (Var) have converged. As expected, 
the mean converges faster than the variance. 

The marginal distributions are approximated by histograms and Gaussian 
fits in Fig. 5. The distributions contain the values of the previous fits of [32] 
and [20] (Tab. 2). 
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Figure 2: Estimated temperature profile of the T40 flyby (left panel) and den¬ 
sities of N 2 and CH 4 for this flyby (right panel). The uncertainties are given at 
a 1 — a level. 



Figure 3: Left: first 10 4 output sample draws for parameters; right: 
(CH 4 ,N 2 ) bimolecular diffusivity calibration. The starting point is determined 
by a maximization of the likelihood function. The rank correlation is —0.85. 
Black cross: [32] fitted values, black plus: [20] fitted values. 


Model 

D(°\ cm 2 s 4 ) 

/? 

Wakeham fit 

0.2019 

1.76 

Massman fit 

0.1892 

1.81 

This work 

0.197 ±0.006 

1.75 ±0.05 


Table 2: Comparison of the parameters’ value. The corresponding Massman 
parameters are given for the Wakeham value using the conversion equation 12 
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Sample index ,^g5 Sample index ,^g5 

Figure 4: Means and variances on the generated samples. 



Figure 5: Posterior histograms. Vertical dashed line: [32] fitted values, vertical 
solid line: [20] fitted values. 
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Figure 6: Estimated bimolecular dif- 
fusivity of (CH 4 ,N 2 ) at 1 atm. The 
dashed lines show the 3 — a interval. 


•nr 2 



Temperature (K) 

Figure 7: Relative uncertainty of the 
bimolecular diffusivity of the couple 
(CH 4 ,N 2 ) at 1 atm. Solid: estimated, 
dashed: measurements. 


4.4 Propagation and convergence of propagated quanti¬ 
ties 

The propagated values of Pch 4 ,n 2 are given in Fig. 6, and its relative uncertainty 
is shown in Fig. 7. The mixture diffusion coefficient of N 2 and its uncertainty 
are given in Fig. 8 and Fig. 9. Fig. 10 shows the convergence of the relative 
uncertainty of the molecular diffusion coefficient of N 2 , which was found to be 
the most difficult to converge, requiring as many as 2 x 10 6 draws. 


5 Discussion 

We have calibrated the Massman parameters for the bimolecular diffusion coef¬ 
ficient for the (CH 4 ,N 2 ) pair using an uninformed prior. The posterior distribu¬ 
tion is given by the product of the likelihood and the prior (Eq. 9). Therefore 
the posterior is informed by data through the likelihood and previously known 
knowledge through the prior distribution. The choice of the prior ensures a 
maximal effect of the information of the data in the Shannon entropy sense. 

As figures 3 and 5 show, the previous fits from [32] and [20] are but one point 
in the posterior space. The cross/dashed lines are [32] values, the plus/solid lines 
are [20] values. Although the fits accurately reproduce the diffusivity at Titan’s 
temperatures, all information about the uncertainty is lost. 

We observe a negative correlation between the relative uncertainty of the 
molecular diffusion of N 2 and the temperature (Fig. 2 and 9). A sensitivity 
analysis of the Massman model (Eq. 13) yields a ln(T/T°) term with respect 
to the j3 parameter that can explain this behavior. Titan’s temperatures being 
lower than the reference temperature of the Massman model (T° = 273.15 K) 
means this log-term will diverge to — oo, therefore increasing the uncertainty 
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Figure 8: Wilke diffusivity of N 2 , the 
3 — a interval is shown with the dash 
lines. 
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Figure 9: Wilke diffusivity relative un¬ 
certainty of N 2 . 
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Figure 10: Convergence of the relative uncertainty of the molecular diffusion 
coefficient of N 2 . 
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of the bimolecular diffusion coefficient (Fig. 7) with respect to /3. A sensitivity 
of the molecular diffusion coefficient with respect to the /? parameter could 
thus explain the negative correlation between the molecular diffusion relative 
uncertainty and the temperature. 


dD a , m = p° fry 

&D(°) P \T°) ’ 

dD s , m _ D (o)P° ( T \ 
d/3 P \T°) 


0 

In 



(13) 


One important consequence is that at lower temperatures, investigations to 
reduce uncertainties on molecular diffusion should focus on the f3 parameter, 
typically by characterizing its value at these temperatures. 

In general, the uncertainties in molecular diffusion are small compared to 
other sources of uncertainty, which can be orders of magnitude larger [13, 4]. 
However, a propagation of the uncertainties in binary molecular diffusion to the 
calculation of the methane mixing ratio as a function of altitude, illustrated in 
Fig. 11, demonstrates that the uncertainty is nontrivial above 1200 km. This 
methane mixing ratio profile is determined using an eddy coefficient that was 
optimized for the methane profile and temperature and density profiles based 
on the T40 flyby as was done in [19]. Also shown are the methane mixing ratios 
measured by the Cassini Ion Neutral Mass Spectrometer [31] during the T40 
flyby. Uncertainties in molecular diffusion coefficients have the greatest impact 
at the highest altitudes, because it is at these altitudes where molecular dif¬ 
fusion has a dominant effect compared to eddy diffusion. This has important 
implications for studies of Titan’s atmosphere that use models to evaluate es¬ 
cape rates based on mixing ratios modeled above 1200 km [7], suggesting that 
uncertainties in molecular diffusion need to be taken into account when making 
conclusions about measured altitude profiles. 

For the common situation where no measurement data is available one has 
to resort to other theoretical or model calculations. In the case of missing labo¬ 
ratory gas mixture data, diffusion coefficients may be estimated by a transport 
theory calculation. The Bayesian framework can also handle this situation; the 
likelihood distribution will, instead of acting on lab measurement data, act on 
transport theory calculation output. It should thus capture the uncertainties 
and inadequacies associated to the theoretical or model calculations. This re¬ 
quires a careful analysis of the corresponding model that is beyond the scope of 
this work. 

In the situations where transport theory calculations yield poor diffusion 
coefficient estimates, it is likely that the uncertainty in the diffusion coefficients 
would be inflated. Access to lab measurements offers the possibility of providing 
favourable uncertainty calculations but investigating the effect of theoretical 
calculations on diffusion coefRent uncertainty is a topic of future work. 
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Figure 11: Propagation of the uncertainties in molecular diffusion, roughly 
±7.5%, to uncertainties in modeled methane mixing ratio. The data points are 
the measurements derived from CASSINI; note that the statistical uncertainties 
from the INMS are orders of magnitude smaller. 
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